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Abstract 

We address the locality problem arising in simulations, which take the 
square root of the staggered fermion determinant as a Boltzmann weight to 
reduce the number of dynamical quark tastes. A definition of such a theory 
necessitates an underlying local fermion operator with the same determinant 
and the corresponding Green's functions to establish causality and unitarity. 
We illustrate this point by studying analytically and numerically the square 
root of the staggered fermion operator. Although it has the correct weight, 
this operator is non-local in the continuum limit. Our work serves as a 
warning that fundamental properties of field theories might be violated 
when employing blindly the square root trick. The question, whether a local 
operator reproducing the square root of the staggered fermion determinant 
exists, is left open. 
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1 Introduction 



Traditionally two formulations of lattice QCD, Wilson [1] and staggered or Kogut- 
Susskind fermions [2,3], were used extensively in practical simulations. Wilson 
type fermions suffer from the breaking of chirality. Staggered fermions rely on 
the appealing idea to get rid of explicit spin degrees of freedom. They arise 
naturally constructing the square root of the lattice Laplace operator [4] . Massless 
staggered fermions have a continuous taste non-singlet U(l) axial symmetry which 
is dynamically broken yielding one Goldstone pion [5-9]. This symmetry prevents 
from additive mass renormalization and together with the reduced number of 
degrees of freedom per lattice site, these properties make staggered fermions well 
suitable for numerical simulations. Their "unwanted" properties are the doubling, 
in the continuum limit each flavor comes in four "tastes", the breaking of the 
taste symmetry at finite lattice spacing and the not straightforward construction 
of operators. 

In addition to Wilson and staggered fermions (and their improved versions) 
many new discretizations of fermions on the lattice have been proposed, most 
notably overlap fermions [10], see Ref. [11] for a recent review. 

With staggered fermions many computations within the quenched approxi- 
mation have been performed [12,13]. An interesting and important universality 
check is represented by the APE plot shown in Ref. [14] where quenched data from 
different fermion discretizations are compared to the continuum extrapolation of 
the CP-PACS collaboration data [15]. 

Since the early times of lattice computations staggered fermions have been 
simulated dynamically [16]. The doubling in four degenerate tastes in the contin- 
uum limit can be exploited since the tastes behave like quarks as far as the QCD 
interaction is concerned. But to describe real world QCD one needs a formulation 
for one or two tastes. An early attempt to describe two tastes by reducing the 
degrees of freedom, called the reduced staggered fermion formalism [17], although 
it satisfies unitarity and positivity, leads to a complex determinant [18] and it 
is therefore not suitable for numerical simulations. Another way of reducing the 
number of tastes is the so called square (or fourth) root trick [19]. It amounts of 
taking the square (or fourth) root of the staggered fermion determinant, motivated 
by its factorization in the continuum limit. While in several places warnings about 
its potential danger have been expressed [19,20], the central problem that this is 
not an ab initio formulation of lattice QCD [11,21,22] has not been addressed in 
the literature. Arguments based on partially quenched chiral perturbation the- 
ory [23] support the square root trick but this should only be considered as a 
motivation for a first principle study. 



2 



From the numerical simulation side there are no major obstacles to implement 
the square root trick. But as more and more results based on such simulations 
are being published and ambitious high precision tests of the standard model 
announced [20], the theoretical issues postponed so far have to be addressed. 
Here we only try to state the problem concerning locality, following a benchmark 
investigation for the overlap case in Ref. [24]. 

At a glance the problem discussed in this work can be introduced by the 
following simple arguments. 

1. In the QCD partition function on the lattice only local operators appear 

Z = [ expi -S g {U) +a A S^i,{x)D^{x) I, (1.1) 

J V&* { x J 

where S g (U) is the gauge action in terms of gauge links U and D is the, 
assumed to be local, Dirac operator acting on fermion fields r/>, if). 

2. The integration over the fermion fields ip, ip generates an effective action 

Z = I det(aD)exp{-S g (U)} = [ exp{-S g {U) +trln(aD)} ,(1.2) 
Ju Ju 

which is non-local in the gauge link variables. 

The non-locality of the effective action does not mean a non-locality of the theory 
if the latter possesses an underlying local formulation in terms of fundamental 
degrees of freedom. 

If the starting point of a lattice theory is eq. (jl.2|) . e.g. with the effective 
action 

S cS = -Sg(U) + hr\n(aM) (1.3) 

corresponding to the Boltzmann weight \J det (aM), then an equivalent formula- 
tion like in eq. (jl.l|) in terms of a local operator D is needed in order to establish 
causality and discuss renormalizability and universality. 

The present work is organized as follows. In Section 2 we formulate precisely 
the locality problem related to the square root trick. We present our "candidate" 
for a local operator D namely the most naive choice, which is to take the square 
root of the staggered operator. In Section 3 we present our analytical results, 
in particular in the free field theory where the square root operator is proven 
to be non-local. Numerical investigations are still of interest to establish the 
actual "non-localization" range in the interacting theory. Section 4 presents the 
simulation results in the quenched approximation. We discuss these results in 
Section 5. There are three appendices devoted to the derivation of our analytical 
results and to the study of finite size effects in our simulations. 
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2 Formulation of the problem 

We consider an Euclidean hypercubic lattice with lattice spacing a and coordinates 
x n — n [i a with n M = 0, 1, . . . , L/a — 1. The number of sites in each direction is 
L/a and it is assumed to be even. The directions are denoted by /i = 0, 1,2,3 
and is a displacement vector by one lattice spacing along the \i direction. The 
action for Kogut-Susskind staggered fermions is 

S = S g (U)-a 4 J2x(x)M X (x) (2.1) 

X 

Mx{x) = mx(x) + ^2^r](x,n)[U{x,fi)x{x + a ll )-W{x-a ll ,fi)x{x-a fl )] , 

where i](x, /i) = (— l)^ v <^ nv are the staggered phases. The fermion field x nas only 
SU(3) color indices. Periodic boundary conditions are used for the gauge field and 
(anti-)periodic for the fermion field. In the continuum limit M exhibits spectrum 
doubling since it describes four degenerate "tastes" , which can be interpreted as 
four flavors of quarks [7,8,17]. At finite lattice spacing there are 0(a 2 ) taste- 
changing interactions even in the free theory. The main advantage of staggered 
fermions is that the action eq. (|2.1j) is invariant for m = under the continuous 
taste non-singlet U(l) axial transformation 

x{x) ^ e **«x(s) , X(x) - e^X(x) , (2.2) 

where 

e( x ) = (_l)E M nM 9 (2.3) 

from which it follows that there is no additive mass renormalization [8,17]. 

The question, whether there exists a formulation for two degenerate tastes 
of staggered fermions has been considered in Ref. [17]. In this so called reduced 
staggered fermion formalism the fermion field \ lives on odd sites only and the 
anti- fermion field \ on even sites only ( "even" and "odd" being defined by the sign 
of e{x) eq. (|2.3|0 . Although this theory satisfies unitarity and positivity it leads 
to a complex fermion determinant [18] and is therefore not suitable for numerical 
Monte Carlo studies. Also this reduced formalism does not have any continuous 
U(l) axial symmetry [8,18]. 

The authors of Ref. [19] proposed to take the square (or fourth) root of the 
fermion determinant to reduce the number of tastes from four to two (or one). 
Based on the factorization in the naive continuum limit 

det (aM) ^4 det(afi) 4 (2.4) 
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in terms of a one-flavor fermion operator fi, taking the square root amounts to 
quenching two out of four tastes and this idea seems to work out in partially 
quenched chiral perturbation theory [23]. 

The square (and fourth) root trick is employed in large scale simulations by 
the MILC [25] and JLQCD [26] collaborations. From the algorithmic side the 
square root of the determinant can be dealt with the R algorithm [16], with the 
PHMC algorithm [26,27] or with the RHMC algorithm [28]. In the computation 
of quark propagators the original four taste staggered operator is used. Correla- 
tion functions are built using sources that project onto the desired valence taste 
components [21]. This is justified in partially quenched chiral perturbation the- 
ory with some complication for taste-singlet operators [23]. However the danger 
with such an approach are possible unitarity violations through a mismatch of 
sea and valence quark lines. This fundamental question, which deserves more 
investigation, is not addressed in the present work. 

In order to establish causality and universality for the theory implicitly de- 
fined by taking the square root of the staggered fermion determinant, a local 
definition of the theory from first principles like in eq. (jl.lj) is needed. Explicitly 
the problem is to find a fermion operator D such that 

det(aD) = Vdet (aM) and \\G(x, y) \\ < Ce^ x - y ^ la , (2.5) 

with C and 7 > independent of U [24,29]. In eq. (|2.5|) y)\\ is the operator 

norm of the kernel 

aDiP(x) = a A Y,G{x,y)il>{y) (2.6) 
y 

and || x — y || e is the Euclidean norm. 

A technical point is that so far we have used M as the operator describing 
four tastes. In numerical simulations the fermion determinant is represented by 
pseudofermions and for this the Hermitean positive definite operator M^M is 
needed, which describes eight tastes. The number of tastes can be reduced to the 
original four by noting that M^M decouples the even and the odd sublattices and 
that [30] 

det(aM) = det(a 2 (M t M) e ) = det(a 2 (M t M) ) . (2.7) 

In eq. (H2J) the subscripts for (M f M) e and (M f M) refer to the operator M+M 
acting on fields living only on the even or odd sublattices respectively. In the 
following to describe four tastes of staggered fermions the operator (M^M) e is 
used. 
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In this work we investigate the locality properties of the so far only known 
candidate for an operator D to satisfy eq. (|2.5jl . namely 

D = y/(MW) e , (2.8) 

where the operator square root is obtained by a Chebyshev polynomial approxi- 
mation [31]. It approximates the unique [32] Hermitean positive definite square 
root of (M+Af) e . 



3 Analytical results 
3.1 Bound 

To derive a bound on the locality of the operator D defined in eq. f!2.8j) we need 
to know the spectral bounds of M^M [33]. In the free theory 

a 2 M^Me ipx = [{am) 2 + sm 2 { Pfl a)}e ipx (3.1) 

which implies for the spectrum spec(a 2 M^M) C [(am) 2 , 4 + (am) 2 ]. In the inter- 
acting case the upper bound is lifted 1 



spec 



(a 2 M f M) C [u,v], u= {am) 2 , v = 16 + (am) 2 . (3.2) 



The operator {M' f M) e has the same spectral bounds [30]. 

In order to handle analytically the operator D we would like to use the series 
of polynomials S n {z) of degree n defined by the generating function 

oo 

Vl+t 2 - 2tz = ^S n {z)t n . (3.3) 

n=0 

Properties of this series (where z is a number) are discussed in the Appendix 
lAl Following Ref. [24] we define the operator (with obvious insertions of unit 
matrices) 

v + u-2a 2 {M^M) e 
v — u 

which maps the spectrum [u,v] of a 2 {M^M) e into [—1, 1]. We then set 

t = e~ e with cosh# = ^^, 6>0. (3.5) 

V — u 



lr The upper bound is equivalent to the one for naive fermions. 
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For this choice of t we have 

4 e -0 

l+t 2 -2tz = a 2 (M f M) e (3.6) 

v — u 

and we can use eq. ()3.3)1 to provide a series expansion for the operator D. 

The kernel G(x,y) eq. ()2.6|) of the square root operator can be expressed in 
terms of the kernels G n (x,y) for the polynomial operators S n (z). From eq. ()3.3|) 
and eq. (|3.5|) we get 

i 00 

G(x,y) = _ v ^Ie e / 2 ^G n (x,y)e-^ (3.7) 

n=0 

Since z as it is defined in eq. (|3.4j) connects two neighboring even sites on the 
lattice, we have 

G n (x, y) = for n < n min = — — , 

2a 

where ||x — y||i/a is the number of links between x and y ( "taxi driver distance"). 
In eq. (jA.9|) a bound for the remainder of the polynomial series in eq. ()3.3|) trun- 
cated after n = n min — 1 terms is given. Since the spectrum of z eq. ()3.4j) is 
contained in [—1, 1] and using the operator calculus, the bound eq. (jA.9|) implies 



a 4 \\G(x,y)\\ < e e/2 e -\\ x - yW a/e) nmin > 2 (38) 

("-min - 1J7T 

for the norm of the kernel G(x,y) in SU(3) color space. Inserting the values for 
the spectral bounds given in eq. ()3.2|) into the definition of 9 eq. ()3.5|) we obtain 

(1771 

9 = ^— + 0((am) 3 ) . (3.9) 

This means that a 4 ||G(x, y)\ is bounded by an exponential oc exp(— r/r loc ) with 
r = ||x — and 

rioc = (3.10) 

The localization range r\ oc stays finite in the continuum limit. This is in contra- 
diction to a local theory where we need r loc proportional to the lattice spacing. 

The result eq. ()3.8|) represents an upper bound on the localization of the 
operator D, the theory could in principle still be local in the continuum. We 
devote therefore our attention to the free theory where we can achieve an exact 
result. Since in the free theory 9 = am + 0((am) 3 ), the bound eq. (J3.8|) gives 
r ioc e) - 2/m. 
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3.2 Exact result in the free theory 

We consider an infinite lattice. The free operator M^M decouples on sixteen 
sublattices of lattice spacing 2a where it acts like a Laplace operator. The Fourier 
transformation of the square root of the diagonal operator in momentum space 
eq. (j3.1J) yields for the kernel G(x,y) eq. ()2.6|) 



where (x M — y^)/a is even for all n, that is x and y live on one of the sixteen 
sublattices of lattice spacing 2a. Clearly then eq. (|3.11|) applies also for D in 
eq. ([2.8)1 . This integral is solved in Appendix |B] and the result (for y — 0) is given 
by replacing a with 2a in eq. ()B.3|) . 

The continuum version of eq. (|3.11|) is also computed in Appendix [0 the 
result is eq. ()B.8|) 



As is shown in Appendix El the continuum and lattice results agree at large Eu- 
clidean distance and this establishes that the operator D in eq. (|2.8J) is non-local 
in the free continuum limit. The localization range is r\ oc = 1/m. By noting 
that ||x||e < ||x||i < 2||x||e we see that the bound eq. (|3.8|) in the free theory is 
saturated. 

It is very unlikely that the introduction of gauge interactions changes qual- 
itatively this result. From the numerical point of view it is still interesting to 
investigate what actually happens in the interacting theory. Because of confine- 
ment the quark mass will be "replaced" in the bound eq. (|3.1U|) by some hadronic 
mass, which could be very large in a favorable case. From a similar study in the 
overlap case [24] it is known that the analytical bound is only poorly saturated. 

4 Numerical results 

To test numerically the locality of the operator D eq. ()2.8|) we follow closely 
Ref. [24]. We consider hypercubic lattices with L/a sites in each direction and 
periodic boundary conditions. The operator D acts on fields living on the even 




(3.11) 




(3.12) 
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/3 


am 


r m G 


L/a 


e 




pol. degree 


#configs. 


6.0 


0.01 


1.29 [34] 


16 


2 x 10" 


-7 


490 


200 








24,32 


8 x 10" 


-9 


735 


600 


6.2 


0.007 


1.34(4) [35] 


24 


4 x 10- 


-9 


1050 


150 








36 


4 x 10- 


-9 


1050 


300 


6.5 


0.005 


1.27(1) [36] 


32 


2 x 10- 


-9 


1470 


300 








48 


2 x 10- 


-9 


1470 


192 



Table 1: Simulation parameters. The last three columns tabulate the relative 
accuracy e eq. (|4.5jl , the degree of the polynomial and the number of configurations 
generated. 

sites only. We define the source field 

= J 1 Hx = yandc = l 
c \ otherwise , 

where y is the location of the source and c runs over the color index of the field. 
We take a point-like color source because we are only interested in the decay 
properties of 

ij){x) = aDU%) (4-2) 

described in terms of the function 

f(r) = max.{\\if)(x)\\\\\x-y\\ 1 = r}. (4.3) 

A different choice of the source in eq. (j4.1j) . e.g. spread in a 2 4 hypercube describing 
"one" taste 2 will not change the result, which concerns a mathematical property 
of the operator that needs to be fulfilled. In eq. f)4.3|) HV'fV)!! is the SU(3) color 
norm and the "taxi driver distance" is defined as 

[|ac — = ^min-fla^ - y^L - \x^ - y^} . (4.4) 

Hence y\\i/a is the number of links for the shortest path on the lattice between 
x and y using the periodicity. The maximal value that — y\\± can take is 2L. 
The taxi driver distance is used here since it arises naturally in the derivation of 
the analytic bound eq. ()3.10|) . 

We compute the function f(r) on configurations generated in the quenched 
approximation and denote the average in the quenched ensemble by (f(r)). We 

2 At finite lattice spacing there are taste changing interactions. 
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Figure V. Decay of (f(r)) defined in eq. f!4.3|) . The curves are normalized to be 
1 at the distance r re f/r = 3.72 corresponding to the location of the rightmost 
vertical dotted line. 



keep the location of the source fixed at y = 0, the average over gauge configu- 
rations projects onto the gauge invariant part of eq. ()4.3|) . Table ^ summarizes 
the simulation parameters. We simulate at three (3 values in order to take the 
continuum limit. The quark masses are obtained from the literature [34-36] and 
define a line of constant physics where the mass tuq of the Goldstone pion ttq is 
roiTiQ = 1.30(3) in units of the hadronic scale r [37,38]. We have two approx- 
imately matched physical volumes Lttlq « 4 and Luiq « 6 at all (3 values. At 
(3 = 6.0 we also have a third larger volume, which we used to study the finite 
volume effects. Statistical errors of derived quantities were determined by the 
method of Ref. [39]. 

We implemented the Chebyshev polynomial approximation using the Clen- 
shaw's recurrence formula [31]. As can be seen from Table ^ for most of the 
computations a relative accuracy e = 10 -8 — 10~ 9 was required, which is defined 
as 

\\(D 2 — (M^M) e )R\\ 



\\(MW) e R\\ 

where R is a normalized Gaussian random vector R. 



(4.5) 
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□ L/a=1 6 
o L/a=24 
x L/a=32 



70 



Figure 2: Effective iocaiization ranges r\ oc (r) at (3 = 6.0, am = 0.01 for different 
volumes. The vertical dotted lines correspond to distances from the source smaller 
than the minimal distances r min (L) L at which finite volume effects for 
L/a = 16, 24 become sizeable. 



In Fig. Qwe show in a semilogarithmic plot the results for (/(r)) as function 
of the taxi driver distance r in units of ro at the three different (3 values for 
the volume LrriQ ~ 6. For a better comparison of the curves we normalized 
(f(r)) by the (linearly interpolated) value (/(r re f)), where the physical distance 
r rof/ r o — 3.72 corresponds to the rightmost vertical dotted line. Fig. [T] shows 
remarkable scaling as the continuum limit is approached, which means that there 
is a physical scale of "non-localization" r\ oc . 

For small distances r <C Aq^ d the form of (/(r)) is dictated by perturbation 
theory and it is expected to follow the polynomial (in 1/r) result in the free theory 
eq. (|3.12|) . For large distances r we assume 



(f(r)) 



oc e 



-r/n oc {r) 



(4.6) 



The inverse localization range ri oc (r) _1 can be computed by taking two consecutive 
values of r to extract the exponential. The results for three volumes at /3 = 6.0 
are shown in Fig- El It is clear that the decay of (f(r)) is not described by a single 
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1 - 
0.5 

I 1 1 1 

0.05 0.1 0.15 0.2 

Figure 3: The continuum limit of the upper bound on the localization range in 
the volume Ltuq ~ 6. 

exponential. On each of the volumes r\ oc (r) has a "bump" for large taxi driver 
distances, which becomes higher as the volume gets larger. This is a finite volume 
effect and is discussed in Appendix IU1 

We adopt two strategies to define a localization range of the operator D in 
the continuum limit. The first is to take for the physical volume LrriQ fa 6 the 
largest value of ri oc (r), which we denote by r™ x , at each j3 value. Changing the 
volume will slightly change this value but the volume LrrtQ « 6 can be considered 
typical for present quenched simulations. The results are shown in Fig. Eland give 
a continuum limit value 

„max 

= 2.15(42) for Lm G «6, (4.7) 

or r™ x m G = 2.8(6). Compared to the analytic bound given in eq. ()3.10|) the 
values of rj^ that we obtain at the three (3 values are about 40 times smaller. 
The fact that the analytic bound is poorly saturated is in agreement with similar 
observations made in Ref. [24]. 

The second strategy is to compute the localization range at a physical value 
of the taxi driver distance from the source. We take two distances, r/ro = 1.86 
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Figure 4: The continuum limit of the localization range computed at two constant 
physical distances from the source. 



and twice this value r/r = 3.72, which are marked in Fig. [TJ Fig. |2]and Fig. [7| 

by the vertical dotted lines. At (3 = 6.0 there are no finite size effects in r\ oc (r) 
computed at these two distances on lattices of size L/a = 24 (Lmc ~ 6), as can 
be seen by comparing to results on the larger volume L/a = 32 in Fig. [7| For 
r ioc( r / r o — 1-86) we checked at the other f3 values that already in the volume 
LvfiQ ~ 4 there are no finite size effects. We take the lattice volume LrriQ « 6 and 
perform a linear interpolation of n oc (r) to get its values at the desired distances 
for the three (3 values. We obtain the results shown in Fig.|3]giving the continuum 
limit values 



n oc (r/r = 1.86) 



r 



n oc {r/r = 3.72) 



r 



0.350(5) , 
0.819(35) 



(4.8) 
(4.9) 



or r\ oc (r/r Q = 1.86)mc = 0.455(12) and r\ 
remind that a local operator would imply r\ oc 



(r/r = 3.72)m G = 1.06(5). We 
r) = 0(a) for all distances r. 
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5 Conclusions 

We have studied the locality problem for a theory defined by taking the square 
root of the staggered fermion operator. In terms of the fundamental fermion fields, 
this operator needs to be local in order to prove causality and unitarity of the so 
defined theory In our work, we proved that in the free field limit such a theory is 
non-local at the scale of the inverse quark mass. Adding gauge fields, simulations 
in the quenched approximation have revealed that the "non-localization" range is 
of the order of the inverse Goldstone pion mass in the continuum limit. We do 
not expect that taking the square root of staggered fermion operators constructed 
with improved actions such as Asqtad [25] or HYP [40] will change qualitatively 
and even quantitatively the situation, since smearing makes configurations even 
closer to the free case. 

How can we interpret our results, also in the context of present day simula- 
tions? We have studied a candidate theory for two tastes of staggered fermions 
obtained by taking the square root of the staggered fermion operator in the ac- 
tion and the corresponding Green's functions. Our work shows that such a theory 
is non-local in the continuum limit at the scale of the Goldstone pion Compton 
wavelength, i.e. the lightest particle in the theory. Hence this solution of the "two 
taste problem" leads to an unacceptable field theory that has to be dismissed. 
This leaves the question open, of course, whether there exists a local operator 
D, which provides a Boltzmann weight equal to the square root of the staggered 
fermion determinant. 

Present dynamical staggered fermion simulations do not use this setup since 
only the square root of the determinant is employed. The corresponding action 
S = l/2tr ln(M 1 'M) e is simulated with the R algorithm [16] with the zero step- 
size limit or S = — 0^((M^M) e ) _1//2 in terms of a psewtMermion field with an 
exact algorithm [26,28]. We emphasize that these simulations generate the correct 
Boltzmann weight since the non-local bosonic formulation is only a technical trick 
to perform the numerical simulations. If a local operator D reproducing the 
square root of the staggered fermion determinant is found then the configurations 
generated by present algorithms are safe. 

There still remains, however, the problem of unitarity. The point is a mis- 
match between the operators used for valence and sea quarks. The valence quarks 
are discretized through the four taste staggered operator, while the sea quarks 
would be introduced through the "still to be found" two taste local operator D. 
We believe that at this point the recovery of the optical theorem is an open prob- 
lem, which needs a clarification. The operator D, acting on the fundamental 
fermion fields, should dictate which are the appropriate Green's functions for a 
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two taste theory. The Green's functions that are used at present are most likely 
not the ones that correspond to the local operator D. 
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A A polynomial expansion of the square root 

Consider the generating function of the Gegenbauer polynomials C%(z) 

oo 

{i + e-2tz)^ = jr,rq>(z) (Ai) 

n=0 

They are defined for any (fixed) 7 G C. Their natural domain is z G [—1,1]. Many 
details about Gegenbauer polynomials can be found in [41], section 10.9. We only 
list the recursion and a bound: 

(n+l)Q +1 W + (n + 2 7 - lJCiUz) = 2(n + 7) zCP n (z) (A.2) 
ICJWI < C^ 2 -'- 1 ) 7 >0 (A.3) 



n 

The power series (|A.1|) converges for \t\ < 1. This may be used as a starting 
point for fractional inversion [42]. Here we specialize to the case of 7 = —1/2, 
which provides an expansion of the square root 



Vl + t 2 - 2tz = ^t n S n {z) (A.4) 

n=0 

Sn(z) = C- l '\ Z ) 

The first few polynomials S n (z) are 

S (z) = 1 
Si(z) = —z 

We find ^(±1) = for n > 2. In fact the recursion (jA.2|) allows us to prove the 
useful relation 



S n (z) = ±-^C^ 3 (z) n>2 (A.5) 
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As an immediate consequence, the bound (|A.3|) provides the estimate 

\S n (z)\ < l - n>2 (A.6) 



In applications of the expansion (|A.4|) . there is particular interest in the 
remainder Q n (z) of the series truncated after the term ~ t n 

oo 

n n ( z ) = tkSk ^ ( A - 7 ) 

k = n + l 

A simple (uniform) estimate follows immediately with the aid of (|A.6|) : 

i 00 1 \t\ n+1 

Pn{z)\ < 2 E 1^=2 1^ (A ' 8) 
k=n+l ' ' 

The coefficient of |£| n+1 diverges as \t\ — > 1. 

We will prove the stronger (and more realistic) bound 

\Q n (z)\ < ^- |t| n+1 n>l (A.9) 

for ze [-1,1] andt G (-1,1). 

Assume t G (0, 1) and let z = cos (p. The relation (|A.5|) translates the integral 
representation of C„ ( [41], section 10.9, eq.(31)) into 

S n (cos(f)) = sin 2 4> / — sin 2 <p (cos</> + i sin 0cos<p) n_2 n>2 
Jo T 

This allows us to perform the summation in Q n (z) eq. ()A.7)) (if n > 1) 

^ , ,. . 9 , „,, H din . 9 (cos</> + 2 sin cosy)™- 1 

Q n (cos0) = sin 2 0r +1 / sin V v 



7r 1 — t(cos + z sin <p cos 99) 



At the end points of the integration, the denominator takes values 1 — te ±l( ^ = 
pe TlS , i.e. we reparametrize (t, <fi) — > (p, 5), with p > and 5 G [0, 7r/2) 



p = |l - £ e ^| = y/\ + t 2 - 2tcos(f) 

cos 5 = (l-tcos0)/p (A.10) 

sin 5 = (t sin 0) /p 

! . 2 d(p . 2 (cos0 + z sin cosy)"' 1 

i2 n (cos0) = t psm — sin ip 

Jo 7T cos — % sin cos <p 
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Rewrite this as a contour integral over u: 



u = 



=>- 1 — pu 
fLfcos d>) 



cos S — % sin 5 cos ip 
t(cos + i sin cos </?) 



/ — \f\ + w 2 — 2u cos 5 (1 — p-u) 
J e -i8 mu 

Je- i 



n-l 



f >, s ^-Q(u)(l-pu) 



with Q(u) 



yl + u 2 — 2u cos 5 
y/{u- e i5 )(u - e~ i5 ) 



The cut of the square root is chosen on the negative real axis. Q(u) has branch 
points at u — e ±tS , with vertical cuts going out to ±ioo. The contour of integration 
connects the two branch points across the real axis. 
Integrate by parts, using Q(u) = at the end points: 



Q n (cos (j>) 



du f d Q(u) 



is iri \du u 
This integral is finally evaluated along the circle 

1 + e 2 ^ 



'l-puj 



(A.ll) 



u{ip) 

Q(u) 
Q(u) 



2 cos 5 



COS ip 

cos 5 



if e [-5, 5] 



coso 



cos^ 



cos 2 5 



a/ cos 2 ip — cos 2 5 



u 



COS if 



fl n (cos(f>) = 



nm 



dip 



-8 



d a/ cos 2 ip — cos 2 5 
dip cos ip 



pu{ip)\ 



Along the contour of integration, 



|1 — pu{ip)\ < t 

Furthermore, Q{u)/u is an even function of <p, which vanishes at the end points 
and reaches a maximum of sin 5 in the middle. In this way we estimate 



|f2„(cos0)| < — / 



dip 



d a/ cos 2 ip — cos 2 5 
dip cos ip 



< 



2t n 



nn 



sin 5 



17 



Inspection of (jA.lOJ) shows that sin 5 < t, which proves (|A.9J) for t G (0, 1). 

The case of t G (—1,0) follows due to the invariance of ft n (z) w.r.t. t — > 
— t, z — > —2. 

B The square root in the free case 
B.l Lattice computation 

On an infinite lattice in d dimensions (lattice spacing a), we study functions of 
the operator (matrix) 

m 2 - A (B.l) 

which reads in momentum space: 

(m 2 — A) (p) = m 2 + 2a~ 2 (l — cos ap M ) 

2 , 4 • 2 fl ?V 

= m + > —sin 



a 2 ~~~ 2 
— 7r /a < Pfj, < 7r/a 

Start from the "generalized propagator" with s G C: 

2s r /a eipx 



(m 2 -A)~ s (x) = a 2s f 

J —5 



. 7r/a (27r) d [a 2 m 2 + £ M 2(1- cos ap M )]' 
Use the standard trick (valid for A > and Re s > 0) 



oo 

s-l-tA 



A- s = — dtt s ^e 

r(s) y 

to factorize the p integrations. They lead to modified Bessel functions of the first 
kind In : 



r*/ a dp 

j — -^exp{ip M x M + 2t cos ap^} = a' 1 I Xfi/a (2t) 



=► (m 2 -A)- s (x) = ^y/ o dtt^e-^ 2 ^Hl^ /a (2t) (B.2) 

For i — > oo, I n (2t) ~ t~ 1 / 2 e 2t and the integral converges if m ^ 0. In case of 
m = 0, we need Res < d/2. 
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For t -> 0: I n (2t) ~ t |n| . In case of x 7^ 0, the product of Bessel functions 
vanishes at least linearly in t. As a consequence, eg. ()B.2|) is valid for Res > —1 
and can be evaluated at s = —1/2: 

— d— 1 /"cc 

(x + 0) (B.3) 

(Note: r(-l/2) = -2y/¥). 

For the special case of x = 0, continuation to s = —1/2 requires to subtract 
a term in the integrand: 

I (2t) d = [l (2t) d -l] +1 

dtt^V^ 2 ™ 2 ^ = r(s)(a 2 m 2 + 2rf)- s 

(m 2 -A) 1/2 (0) = a~ d - 1 (a 2 m 2 + 2rf) 1 / 2 (B.4) 

POO 

/ dtr 3 / 2 e-'( a2m2+M ) [7 (2t) d - 1] 
Jo 



a -d-i r°o 



r(-i/2) 



B.2 Continuum calculation 

Compute the "generalized propagator" again: 

"°° d d p e ipx 



(m 2 -A) s (x) = J 



^ (27r) d (m 2 + p 2 ) s 

1 [°° J++S-1 f ipx-(m 2 +p 2 )t 

T(s)J J (2nY 6 



(47T) 



-d/2 



T(s) 



dtt s-l-d/2 e -m 2 t-x 2 /(4,t) 

(B.5) 



This time, we end up with modified Bessel functions of the second kind K v : 

/\\rr\\ \ s - d / 2 roo 



2m 

x 



s-d/2 



2K d/2 _ s (m\\x\ 



2m 

(A n )-d/2 /II || \ s-d/2 

(m 2 -A)-» = 2K d/2 . s {m\\x\\) (B.6) 
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Figure 5: Comparison of the lattice eq. ()B.3j) and continuum eq. (|B.8|I results 



The derivation is valid for i ^ and < Res < d/2, but analytic continuation 
to s = —1/2 is obvious: 



(m 2 - A) 



1/2 



X) 



(47T) 



-d/2 



T(-l/2) \2m 



-(d+l)/2 



2K [d+1)/2 (m\\x\\) (B.7) 



In four dimensions, we make use of the fact that Bessel functions with half-integral 
index are elementary, e.g. 



K B/2 (z) 
(m 2 - A) 1/2 (x) 



7T\V2 _ / 3 3 

2z) e [ 1 + - Z + 7* 
1 m 2 



An 2 \\x\\ 3 
3 



""' '' ' 1 + - + ° 



4tt 2 11 
(d = 4) 



x \ 1-5 e -m||x|| 



m||x|| m 2 IMI 2 



(l + m| |x| | + m 2 | |x| | 2 /3) 



(B.8) 
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B.3 Comparison 

The result in eq. flB.8|) is the (formal) continuum limit of the lattice expression 
eq. (|B.3|) . We expect 

In fact, it has been verified numerically that the lattice expression approaches the 
continuum result with increasing ||x||/a, using the Euclidean distance \\x\\e on 
the lattice as well. This is shown in Fig. El for am = 0.05. 

C Finite volume effects 

When applying the staggered operator M^M hops to neighbors are suppressed 
by a factor l/(4(am) 2 + 8) with respect to the static mass term. So aD£ c (x) 
in eq. (|4.2|) receives smaller contributions from paths with a larger number of 
hops connecting x with the source at y. As the taxi driver distance \\x — y\\\ 
eq. (|4.4j) grows, the relative weight of path wrapping around the lattice grows. 
This growth though is expected to be smaller in comparison on larger lattice sizes 
V > L simply because the path "around the world" is longer. 



y 



"J c 




9 ° y c ' x ° 9, '- 




X 





L 



y x 

Figure 6: Representation in one dimension of a finite volume effect. The thickness 
of the lines representing the hops is proportional to the weight of their contribution 
in eq. (|4.2|) . The vertical dotted lines mark the largest possible taxi driver distance 
from the source at y. 

The situation is schematically represented in Fig. |H1 Thus we expect 

(f(r))\ L >(f(r))\ L ,, L'>L, r>somer min . (C.l) 

This expectation is confirmed by the numerical results shown in Fig. which are 
obtained at (3 = 6.0, am = 0.01 comparing (f(r)) computed on I/a = 24 with 
V ja = 32 lattices for < r < 2L. We observe that the minimal distance at which 
the finite size effects become sizeable is r m j n pa L. The vertical dotted lines in 
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Figure 7: Finite volume effect at (3 = 6.0, am = 0.01. Comparison of (f(r)} 
computed on L/a = 24 with L'/a = 32 lattices for < r < 2L. 



Fig. [7| correspond to the distances used in eq. (|4.8|) and eq. (|4.9|) . At these two 
values of r finite volume effects are absent on the L/a = 24 lattice. 

The effect leading to eq. (jC.l|) is the dominant (but not the only one) finite 
volume effect. It is responsible for the "bumps" of the effective localization range 
at large distances r visible in Fig. EJ which vanish as the volume increases. The 
size of this finite volume effect decreases as the mass am increases. 

Concluding this section we remark that the maximization operation in the 
definition of f(r) eq. (j4.3J) solves the ambiguity for the case when points at the 
same taxi driver distance but different Euclidean distances from the source are 
considered. The maximum is presumably attained at the point with the smallest 
Euclidean distance, which is also the one least affected by the four-dimensional 
generalization of the finite size effect depicted here in Fig. |U1 
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